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The construction of exponentially localized Wannier functions for a set of bands requires a choice 
of Bloch-like functions that span the space of the bands in question, and are smooth and periodic 
functions of k in the entire Brillouin zone. For bands with nontrivial topology, such smooth Bloch 
functions can only be chosen such that they do not respect the symmetries that protect the topology. 

This symmetry breaking is a necessary, but not sufficient condition for smoothness, and, in general, 
finding smooth Bloch functions for topological bands is a complicated task. We present a generic 
technique for finding smooth Bloch functions and constructing exponentially localized Wannier 
functions in the presence of nontrivial topology, given that the net Chern number of the bands in 
question vanishes. The technique is verified against known results in the Kane-Mele model. It is 
then applied to the topological insulator Bi 2 Se 3 , where the topological state is protected by two 
symmetries: time-reversal and inversion. The resultant exponentially localized Wannier functions 
break both these symmetries. Finally, we illustrate how the calculation of the Chern-Simons orbital 
magnetoelectric response is facilitated by the proposed smooth gauge construction. 


I. INTRODUCTION 

Discoveries of the past decade established that not only 
the geometry of electronic bands but also their topology 
is a fundamental material property. Being at the root 
of a variety of physical effects,topologically nontrivial 
bands turned out to be ubiquitous, and many materials 
were found to host various topological states.®^® 

As a general principle, nontrivial topologies in band 
structures are a consequence of some symmetry. They are 
examples of symmetry protected topological states. 

For instance, magnetic insulators in 2D realize the integer 
quantum Hall effect in the absence of external magnetic 
field.jjj such insulators, dubbed Chern insulators, 
the nontrivial topology is captured by a Chern number 
topological invariant,which cannot be changed unless 
the U(l) charge conservation symmetry is broken.^® In 
the case of time-reversal (TR) symmetric insulators, a 
Z 2 topological invariant is assigned^® to the band struc¬ 
ture, that can not be changed without breaking the 
TR-symmetry (provided the band gap remains open). 
This concept of symmetry protection can also be gen¬ 
eralized to crystalline symmetries, which protect the 
value of certain topological invariants, associated with 
these symmetries,giving rise to crystalline topologi¬ 
cal insulators . 20 

The topology of a band structure in all these cases is 
associated with the Bloch bands and their Berry curva- 
ture^^’^^ in momentum space. An alternative description 
of crystalline solids can be made in position space by 
means of localized Wannier functions (WFs).^^’^^ Such a 
local basis is often preferable to that of Bloch functions, 
for example for modeling finite size effects in materials, 
chemical bonding,^® computing electronic polarization,^^ 
or ballistic transport.^® This variety of applications and 
the existence of established numerical techniques for WF- 
based analysis of materials,makes it important to 
extend the phenomenology of WFs to materials with non¬ 


trivial topology. 

Another motivation for finding Wannier representation 
of topological bands is the calculation of the isotropic 
contribution to the orbital magnetoelectric response.®® 
The linear magnetoelectric coupling tensor is defined 
as 


aij 


m\ ^ (dMA 


( 1 ) 


where P and M are the polarization and magnetization of 
a material, E and B are electric and magnetic fields and 
both derivatives are evaluated at zero fields. While lattice 
and spin degrees of freedom also contribute®®’®^ to aij 
here we focus solely on the orbital (frozen-lattice) part. 
The orbital magnetoelectric response can be further split 
into two parts®®’®^ 


^ij T 2,7rh ^ 

of which aij is traceless, and the second, isotropic, part 
is characterized by a dimensionless quantity 9, called 
the “axion angle” in high energy physics.®’®® There are 
two contributions to 0:®^ one is an ordinary perturba¬ 
tive Kubo term and the other is a purely geometrical 
one, the Chern-Simons contribution Ocs- The latter one 
can be evaluated from the ground-state electron wave 
functions, by computing an integral of the Chern-Simons 
3-form over the entire Brillouin zone (BZ).®’®®’®^ How¬ 
ever, this evaluation requires a choice of Bloch states 
that are smooth and periodic in k-space as detailed in 
Sec. VD. Alternative formulations of this term in posi¬ 
tion space require the existence of exponentially localized 
WFs (ELWFs).®® 

Magnetoelectric response was long thought to be ob¬ 
servable only in materials that break TR and inver¬ 
sion symmetries. Remarkably, it turned out that TR- 
and inversion-symmetric topological insulators (TIs) are 
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characterized by a nonzero quantized magnetoelectric 
response. 33.36 presence of one of these symme¬ 

tries the Kubo contribution to d vanishes and Q = 9cs- 
The seeming contradiction is resolved by noticing that 
0 couples to E • B term in the Lagrangian, so that the 
corresponding equations of motion are invariant under 
9^9 + 27r.^® Since, in addition, 9 is odd under both TR 
and inversion,the two values 9 = 0,Tr are compatible 
with these symmetries. 

In TR- or inversion-symmetric TIs these two values of 
0 can be used as an analogue of the Z 2 topological invari¬ 
ant, ^^3 -vvith 9 = TT corresponding to the TI phase and 
0 = 0 to the normal insulator (NI) one. For the cases of 
TR- and inversion-symmetric insulators the quantization 
of 9 is exact and since 9 corresponds to the Z 2 invari¬ 
ant, it can be obtained using the methods for comput¬ 
ing this invariant.However, for materials that lack 
these symmetries, 0 is in general not quantized, and it 
becomes necessary to evaluate the 0cs-term directly. As 
mentioned above, this requires ELWFs. Since nontrivial 
band topologies can arise even when there are no symme¬ 
tries that quantize the magnetoelectric response, compu¬ 
tation of 0CS gives yet another motivation for developing 
a method to obtain ELWFs for topologically nontrivial 
band structures. 

WFs are constructed by Fourier transforming Bloch 
states, and thus the momentum space geometry of the 
Bloch function can strongly influence the properties of 
the resultant WFs, in particular the degree of its local¬ 
ization. An important example is that the construction 
of an ELWF for a Bloch state with nonzero Chern num¬ 
ber is impossible.This is a consequence of the fact that 
a Chern number represents an obstruction for choosing 
the Bloch state to be a smooth and periodic function of 
k globally in the whole BZ.^^ 

ELWFs can only be constructed if all the Bloch states, 
for which a Wannier representation is constructed, are 
smooth and periodic in the whole BZ. For topologically 
trivial bands, the choice of smooth Bloch states, referred 
to as a smooth gauge choice, respecting all the symme¬ 
tries of the underlying band structure is generally possi¬ 
ble. For topological bands, however, the symmetry that 
protects the topology represents an obstruction to choos¬ 
ing a smooth gauge that respects this symmetry. Conse¬ 
quently, the smooth WFs have to to break this symme¬ 
try. This was explicitly illustrated to be the case for the 
time-reversal (TR) symmetric Z 2 

While it is clear that the construction of a smooth 
gauge on a lattice requires breaking certain symmetries 
in the gauge, finding an explicit representation of smooth 
Bloch states on a lattice of k-points, required for a nu¬ 
merical construction of ELWFs, is a nontrivial task. An 
explicit construction based on parallel transport of the 
occupied states was obtained in Ref. 45 for a 2D model 
of a quantum spin Hall insulator. That construction, 
however, is tedious to generalize to the many-band case, 
and especially to higher dimensions. 

Another approach to finding a smooth gauge is based 


on projecting certain localized orbitals that break the 
topology-protecting symmetry onto the the occupied 
states (see Sec. H for details) and is more appropriate 
for material calculations.^^ The problem of this method 
is that no specific algorithm for choosing the orbitals for 
the projection is presented. The requirement of breaking 
the topology-protecting symmetry in the initial projec¬ 
tion is necessary but not sufficient for finding ELWFs. 

In this work we develop a generic algorithm to con¬ 
struct ELWFs for a set of topologically nontrivial bands, 
with a zero net Chern number. The idea is to construct 
an adiabatic connection between trivial and TI phases by 
breaking the symmetries that protect the topology every¬ 
where along the connecting path, except the initial and 
final points. Since the topology-protecting symmetries 
are broken at intermediate steps, it becomes possible to 
find a path connecting these two phases in a parameter 
space such that the insulating gap remains open along 
the connection. 

The construction of this path allows one to avoid the 
problem of finding the correct symmetry-breaking pro¬ 
jection in the topological phase. Instead, one discretizes 
the path in parameter space into steps, and constructs 
ELWFs at each step by projecting onto ELWFs found 
at the previous step. The initial step requires finding 
ELWFs in a topologically trivial case, which is a stan¬ 
dard task.^® As a result, one naturally obtains ELFWs 
for the topological bands at the final point of the dis¬ 
cretized path. 

The breaking of symmetries plays a key role in adiabat- 
ically connecting two topologies. It is often the case that 
there are more than one symmetry protecting the topol¬ 
ogy of the band structure, and it is important that all of 
them are broken along the adiabatic path. For example, 
the Z 2 topology of the TR-symmetric TIs can be addi¬ 
tionally protected by inversion, in-plane mirror^^’^® 
and certain other point group symmetries.^® These sym¬ 
metries need to be broken in addition to TR along the 
adiabatic path to obtain a smooth gauge in such band 
structures. 

This paper is organized as follows. In Sec. H we pro¬ 
vide the theoretical background about construction of 
ELWFs. In this light, we also discuss the topological 
obstruction for constructing ELWFs in TIs. In Sec. HI 
In Sec. IV we apply our technique on the model of Kane 
and Mele^®. Then in Sec. V we apply our technique onto 
Bi 2 Se 3 and use it to calculate the Chern-Simons magne¬ 
toelectric coupling 0CS- Finally, we summarize our find¬ 
ings and give an outlook in Sec VI. 


II. WANNIER FUNCTIONS AND 
TOPOLOGICAL OBSTRUCTION 

A. Construction of Wannier functions 

Here we provide a brief review of the methods intro¬ 
duced in the work of Ref. 26 used for an explicit con- 
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struction of ELWFs starting from a set of Bloch states 
obtained by diagonalizing the Hamiltonian on a mesh 
of k-points. The construction requires a projection of a 
set of trial localized orbitals onto the Bloch states. The 
choice of these orbitals, which is simple in the case of 
topologically trivial bands, becomes problematic for the 
case of nontrivial topology. For the method we present in 
this work, localized trial orbitals need to be found only 
in the topologically trivial case. 

Given a Bloch state ^/>„k(r), a corresponding WF is 
defined as 

(r|Rn) = fF„(r-R) = —^ [ dke"**"'^'!/'«k(r), (3) 

where V is the volume of the unit cell, d is the dimen¬ 
sionality. The Bloch wave functions tjjnk are assumed to 
be normalized within the unit cell. 

This definition, however, is not unique. The non¬ 
uniqueness, referred to as gauge freedom, is easily seen 
when constructing WFs for a set of N Bloch states. A 
general unitary transformation U(A^) of the N occupied 
bands 

It^nk) — E Umn (kjl'^mk), (4) 

ra 

results in a different set of Bloch states that span the 
same Hilbert space as the original ones, and hence can 
equally be used for constructing the Wannier represen¬ 
tation. Depending on the particular gauge choice, the 
resultant WFs and their degree of localization can vary 
a lot. In order to obtain FLWFs, the gauge choice has 
to be smooth, meaning that all N Bloch states used to 
construct the WFs are smooth in the entire BZ, and obey 
periodic boundary conditions ipnk = V'nk+G upon trans¬ 
lation by any reciprocal lattice vector G. 

But even such a smooth gauge choice is not unique, 
since a smooth gauge transformation performed on a set 
of smooth Bloch states will result in a different set of 
smooth Bloch states and different ELWFs. Additional 
constraints can be put on the gauge to reduce the gauge 
freedom. A very common choice of such a constraint is 
the requirement of maximal localization of the resultant 
WFs in position space proposed in the work of Ref. 26. 
This gauge is obtained by minimizing the spread func¬ 
tional 

N 

^ = {r)l) (5) 

n—1 


where {r)n = f \Wn\'^rdr, and H/ and fl stand for the 
gauge- indep endent 


N 

^^i = E 

n—1 

and gauge-dependent parts of the spread 

AT 

^ = E E l(R'W|r|On)p. 

n—1 Rm^^On 


(r^)„ - |(RTO|r|On}|^ 
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The goal of maximal localization is to find a U(iV) trans¬ 
formation that, when applied to some initial set of Bloch 
states according to Eq. (4), minimizes H to produce a 
set of maximally localized WFs. Maximal localization 
can always be performed, once the initial choice of Bloch 
states is smooth. 

The necessity for smoothness of the initial Bloch states 
is most easily seen from the general procedure used to 
construct a set of WFs, as described in Ref. 26. Following 
this procedure, to construct a set of N WFs from N 
isolated (that is, separated by energy gaps from the rest 
of the spectrum) Bloch bands obtained from numerical 
diagonalization of the Hamiltonian, a set of N localized 
trial states jr^) is chosen. For each momentum k all of 
the trial orbitals are projected on the Bloch states to get 
a set of N Bloch-like states^® 


N 

|Tik) = Acln) = |V'nk)('i/’nk|Ti), (8) 

n—1 


which are not orthonormal. 

For the construction of WFs these states need to be 
orthonormalized, and it is the orthonormalization proce¬ 
dure, where the smoothness of the gauge becomes crucial. 
To see this, apply Lowdin orthonormalization procedure 
to the states, which is commonly used to get a set of or¬ 
thonormalized Bloch states, to orthonormalize the states 
lUk) 


i^„k) = y][5(k)-i/2]^„|T^k), (9) 

m 


where 


^^„(k) = (T^k|T„k) (10) 


is the overlap matrix. While IV'nk) are not the eigen¬ 
states of the single particle Hamiltonian, they span the 
same space as the usual Bloch eigenstates, and therefore 
describe the same ground state. For the trial states em¬ 
bodying a reasonable assumption about the character of 
the bands described, the IV'nk) will be smooth functions 
of k. In this case the WFs constructed from them by 
means of Eq. (3) are expected to be exponentially local¬ 
ized, and the degree of localization is further increased 
by doing maximal localization. 

However, the orthonormalization procedure breaks 
down if at some k the determinant of the overlap ma¬ 
trix vanishes, that is det5'(k) = 0. This is guaranteed 
to happen if one runs into a topological obstruction, for 
example, by projecting onto a set of trial states In) that 
respect the topology-protecting symmetry, as shown be¬ 
low. 
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FIG. 1. Possible flows of the hybrid Wannier charge centers in 
a 2D TR-symmetric system. Panel (a): topologically trivial 
insulator. Panels (b) and (c): quantum spin Hall insulator. 
A TR-symmetric gauge is used in panels (a) and (b), while a 
TR-breaking gauge is used in panel (c). 


B. Topological obstruction via hybrid Wannier 
functions 

Hybrid WFs®° (HWFs) differ from usual WFs in that 
the Wannier decomposition is done in one direction only 

\Ryka:n) =-^ [ dkye~"^«^''\'ipnk). ( 11 ) 

J-T7/a 

This wave function is localized in momentum in the x- 
direction and in position in the y-direction, being a hy¬ 
brid of Bloch- and Wannier-like functions. HWFs proved 
to be useful for classifying topologies of band struc¬ 
tures,^®’®^ providing an intuitive approach to topological 
invariants. 

The band structure topology can be obtained by track¬ 
ing the charge centers of HWFs defined as 

Vnikx) = {0kxn\y\0kxn), (12) 

where the HWFs are located within the home unit cell. 
Computation of the charge centers does not require an ex¬ 
plicit construction of HWFs, and can be done by means of 
the parallel transport procedure, as described in Ref. 38. 
The HWF charge centers obtained in this gauge®^ are the 
eigenvalues of the projected position operator.®® How¬ 
ever, for an isolated set of N bands the gauge can be 
chosen differently, resulting in different values of y{kx), 
and it is only the sum of all the N centers that is gauge 
invariant (modulo a lattice vector in the y-direction) at 
each /ca;.®® 

The parallel transport gauge respects the symmetries 
of the system. For example, in the presence of TR- 
symmetry the HWF charge centers come in Kramers 
pairs yn{kx) = ym[—kx) being doubly degenerate at the 
TR-invariant momenta —k* = k* + Gx, where Gx is a 
reciprocal lattice vector in the cc-direction. This is illus¬ 
trated in Fig. 1(a) and (b) for a model with two occupied 
bands. The two centers are degenerate a,t kx = 0 and 
kx = ±TTja. Both centers evolve smoothly in between 
±7r/a. In the case of Fig. 1(a) they return to the orig¬ 
inal value at the BZ boundary and the band structure 


is topologically trivial. In the case of Fig. 1(b) they in¬ 
terchange, which is generally the case in a quantum spin 
Hall insulator.®^ 

This interchange results in a discontinuity of the hy¬ 
brid Wannier charge center lines as a function kx'. the 
two centers are continuous functions of momentum for 
kx G [~7r/a,7r/a], but the periodicity constraint is not 
satisfied, meaning that the center position does not nec¬ 
essarily return to the original value after a 27r/a change in 
momentum. This is equivalent to placing the topology- 
dictated gauge discontinuity on the boundary of the BZ. 
In this gauge choice, the Wannier centers each correspond 
to a particular, individual Chern number, since the shift 
of the Wannier center at the boundary is still an integer 
number of unit cells. Hence, each of the hybrid Wannier 
functions can be understood as a well-defined function, 
in a sense that a Chern number can be assigned to each 
of them, despite the possible presence of degeneracies in 
the energy spectrum. 

For example, the two bands in the Kramers pair are 
degenerate at the TR-symmetric momenta, but the cor¬ 
responding hybrid Wannier centers result in a splitting of 
the pair into two subspaces with well-defined Chern num¬ 
bers. For a trivial insulator, illustrated in Fig. 1(a), the 
two centers correspond to zero Chern numbers. For the 
case of a quantum spin Hall insulator the hybrid Wan¬ 
nier centers shown in Fig. 1(b) indicate that in the cho¬ 
sen gauge the occupied space is split into two subspaces 
with Chern numbers equal ±1.®^ In both cases the two 
subspaces are TR-images of each other, but the nonzero 
Chern numbers of the two states in the topological phase 
signal that the TR-symmetric gauge is not smooth on 
the whole BZ torus. The above argument is equally valid 
for more than two occupied bands, where each of the two 
subspace contains several bands and Chern numbers are 
then assigned to each subspace. 

The Kramers degeneracy of the Wannier centers per¬ 
sists for any gauge, in which TR maps one state at k onto 
the other state at —k.'*® For this reason TR-symmetry in 
the gauge represents an obstruction for it to be smooth 
and periodic in the BZ. Tracking the presence of this 
obstruction is the basis for the methods of computing 
topological invariants. ®^^®® 

The above analysis, in accord with other meth¬ 
ods,leads to the conclusion that in the Z 2 TIs 
a smooth gauge, and hence ELWFs, have to break TR- 
symmetry. Breaking the symmetry in the gauge (but 
not in the Hamiltonian) lifts the Kramers degeneracy 
of Wannier centers at the TR invariant momenta, and 
can result in a smooth gauge^® like the one illustrated in 
Fig. 1(c). The Wannier centers are smooth and single¬ 
valued throughout the BZ. This means that the occupied 
subspace was split ^® into two subspaces with zero Chern 
numbers, and hence the construction of the ELWFs be¬ 
comes possible in this gauge. 

The above line of reasoning can be easily generalized 
to band structures, where the nontrivial topology is pro¬ 
tected by a symmetry different from TR, for instance a 
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Htopological 



FIG. 2. Schematic path in a parameter space of some model 
Hamiltonian. The adiabatic path is chosen to connect the 
topologically trivial and nontrivial phases avoiding gap clo¬ 
sures. The topology-protecting symmetries is broken along 
the path. 

crystalline A particular example can 

be that of the crystalline TI SnTe,®h9 where the mirror 
Chern numbers^^ of ±2 can be defined on the {110} mir¬ 
ror plane in the BZ. Hence, the mirror symmetry should 
be broken in the gauge to obtain ELWFs for this ma¬ 
terial. Again, not any symmetry-breaking gauge would 
work, but only specific ones that provide a decomposi¬ 
tion of the occupied space into states that are smooth in 
the BZ. 


III. AVOIDING TOPOLOGICAL 
OBSTRUCTION 

Finding suitable initial projections, for which 
detS'(k) 0 in the whole BZ, is particularly difficult 
for TIs due to the band inversion and hence, band char¬ 
acter change, present in these materials. According to 
the discussion above and also the previous findings,^® the 
initial projections must be chosen such that they break 
the topology-protecting symmetries. However, as men¬ 
tioned above, this requirement is necessary, but not suffi¬ 
cient. Given that with an increasing number of bands the 
number of possible projections increases exponentially, a 
brute force search for suitable projections is complicated 
and inefficient. 

Here we describe a technique for obtaining a smooth 
gauge and ELWFs for topological bands that avoids the 
need of finding suitable initial projections. Instead, we 
construct an adiabatic path that connects a topologically 
trivial Hamiltonian to the topological one, such that the 
system remains gapped along the path, as illustrated in 
Fig. 2. To keep the system gapped, all the symmetries 
that protect the topology need to be broken along the 
path. The initial topologically trivial Hamiltonian can be 
chosen in different ways. For example, in calculation for 
a material, where the nontrivial topology is often driven 
by spin-orbit coupling (SOC), the Hamiltonian with SOC 
artificially set to zero is topologically trivial and can be 
used as a starting point, if it is gapped. 


Finding suitable trial functions |ri) for topologically 
trivial bands is usually straightforward, since they re¬ 
spect the symmetries of the Hamiltonian.^® Moreover, an 
automated procedure for finding optimal projections for 
the trivial bands was proposed recently in Ref. 56. Thus, 
the ELWFs at the initial step of the path can always be 
found. 

At least two parameters are required to parametrize 
the adiabatic path. One controls the topological phase 
transition by tuning the strength of symmetry breaking 
along the path.®^ Care should be taken when breaking the 
symmetries, since double symmetry protection can occur, 
meaning that there is more than one symmetry protect¬ 
ing the nontrivial topology. The issue of double symme¬ 
try protection is discussed with more detail in HI A. 

Once the adiabatic path is constructed, the smooth 
Bloch states for the initial, topologically trivial, Hamil¬ 
tonian TLqj are found and used to construct ELWFs. Af¬ 
ter this is done the path is discretized into L steps, so 
that Hl is the Hamiltonian of the TI in question. At 
step I the Hamiltonian Hi is diagonalized and the EL¬ 
WFs obtained at the initial step are used as the trial 
orbitals for its occupied states. For dense enough dis¬ 
cretization, the corresponding projection is guaranteed to 
have det5'(k) ^ 0 throughout the BZ, since Hi and Hq 
are only slightly different. The resultant states are used 
to construct ELWFs for Hi. The procedure continues by 
projecting the ELWFs obtained at the step £ onto the 
occupied state of the Hamiltonian until the final 

point of the path is reached. Since at each step ELWFs 
were found, one ends up with ELWF-representation of 
the occupied topologically nontrivial space, which solves 
the problem of finding a smooth gauge.®® 

We emphasize that this method is general, being appli¬ 
cable to any isolated set of bands and in any dimension, 
provided that the net Chern number of this set is zero. 


A. Double Symmetry Protection 

Here we discuss the double symmetry protection of the 
topological phase using two examples: Z 2 TR-symmetric 
insulators with mirror or inversion symmetries. 

We start by considering the case of coexisting mirror 
and TR symmetries. When a TR-symmetric plane in the 
BZ of the insulator is invariant under mirror symmetry, 
the nontrivial 2D Z 2 invariant of this plane suggests that 
both TR and mirror symmetries need to be broken in 
the smooth gauge. The mirror symmetry breaking in the 
gauge means that the Bloch states |V^„k) of Eq. 4 obtained 
in the smooth gauge and used to construct ELWFs are 
not eigenfunctions of the mirror operator on the mirror 
plane. This can be seen by noting that the two states in 
a mirror-symmetric Kramers pair have opposite mirror 
eigenvalues ±1 on this plane, so the occupied space of 
an insulator on this plane can be split into two subspaces 
according to the value of the mirror eigenvalue. In the Z 2 - 
odd phase +i and —i have opposite odd Chern numbers. 
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and thus breaking TR symmetry only is not sufficient to 
remove the topological obstruction, since the two mirror- 
labeled subspaces remain nontrivial. Hence, the mirror 
symmetry has to be broken as well, to construct Bloch 
states that are smooth in the entire BZ for the Z 2 -odd 
phase. 

Another example of double protection is that of an 
inversion- and TR-symmetric TI. The smooth gauge in 
this case also has to break both symmetries, as we ar¬ 
gue below. Several works^®’^® discussed the possibility of 
a topological classification in the presence of inversion- 
symmetry only. In particular, the work of Ref. 46 pre¬ 
sented a study of the HWF centers in inversion symmet¬ 
ric insulators, obtained by a special construction (called 
a Wilson loop^®). In that construction the projector onto 
all the occupied states is defined at each fc-point, and 
the HWF centers are obtained by taking the log of the 
eigenvalues of the product of operators Hk taken over 
discretized values of A: on a closed loop in momentum 
space^®. For inversion symmetric systems the inversion 
symmetry I puts the following constraint on the projec¬ 
tor P° = iP°J. 

The resultant HWFs are inversion symmetric in a sense 
that for each center yn{kx) there exists an inversion sym¬ 
metric partner —yni—kx), where n is not necessarily 
equal to m. Thus, in the presence of inversion symme¬ 
try the sum of all the HWF centers at the inversion- 
invariant momenta = {0,7r/a} can only be 0 or 6/2, 
which are the inversion-symmetric values assuming the 
inversion center coincides with the center of the unit cell. 

However, as discussed above, when constructing EL- 
WFs, one needs a representation of the occupied sub¬ 
space in terms of Bloch states that are smooth and pe¬ 
riodic in k-space, meaning that the projector on each of 
these Bloch states is smooth. Thus, in a smooth gauge 
the net projector onto the occupied states is decomposed 
into a set of projectors P^ = which is 

smooth. However, imposing inversion symmetry con¬ 
straint of the form on each of the individ¬ 

ual projectors, restricts the corresponding HWF centers 
(not their sum, but each of them separately) to take on 
inversion-symmetric values at kx = {0,7r/a}, that is ei¬ 
ther 0 or 6/2. We call a gauge, in which each of the 
projectors respects inversion symmetry in this sense an 
inversion-symmetric gauge. 

In an inversion-symmetric gauge each Bloch state 
IV'nk) respects inversion symmetry e®‘^|'!/'nk) = Ilfpri-k), 
and hence the corresponding WF centers are subject 
to the condition (r)„ = — (r)„ mod R. Consequently, 
the HWF centers fulfill yn{kx) = —yn{—kx) modulo 6 
and are restricted to the values yn{kx) = {0,6/2} at 
kx = {0,7r/a}. Note, that the HWFs are smooth in mo¬ 
mentum in the interior of the BZ, and thus the index 
n refers to the same center on both sides of this equa¬ 
tion. If for fc* = {0,7r/a}, jV'nk) has the same parity at 
ky = 0 and n/b then y„(fc}) = 0, while if the parities are 
opposite yn{k%) = bl2. 

This means that if an inversion-symmetric Bloch state 


has a non-zero Chern number on some 2D BZ (2D cut of 
a 3D BZ), thus being not smooth, it is impossible to make 
it smooth (that is, change the Chern number to zero) if 
the inversion symmetry in the gauge is preserved. This is 
shown in Fig. 1. The illustrated centers are obtained in 
an inversion-symmetric gauge, and by construction are 
smooth in the interior of the BZ kx € [—7r/a, tt/o]. The 
Chern numbers of the centers shown as solid blue and 
dashed red lines in Fig. 1(b) are given by the number 
of unit cells traversed by the HWF center when going 
from one edge of the BZ to the other, and are equal 
to ±1. They cannot be changed to zero by breaking 
TR-symmetry alone, while preserving inversion symme¬ 
try in the gauge in the above sense, since that would 
require moving the center position away from the value 
fixed by inversion symmetry at some high-symmetry mo¬ 
mentum. In the inversion-symmetric gauge, where each 
Bloch state separately is taken to be the eigenstate of 
inversion operator, this argument holds for any number 
of Kramers pairs, since the argument applies for each 
individual HWF center. 

To state it more rigorously, the TR and inversion in¬ 
variant Z 2 topological insulating state in 2D is character¬ 
ized by the total number of negative parity eigenvalues of 
occupied states at the four TR invariant momenta:'*®’®® if 
this number divided by two is odd the system is a topo¬ 
logical insulator.®® Let nj" be the difference in the number 
of negative parity eigenvalues between points (0, 0) and 
(0,7r/6), and the difference between points (tt/o, 0) 
and (tt/o, 7r/6). If n/ and are different (which is the 
case in a TI) then the number of HWF centers for which 
Vnikx) = b/2 is also different at = 0 and kx = tt / o . 
Since the HWF centers are smooth in the BZ interior 
kx G [— 7r/a, tt/o], at least |n/ —n^l of the hybrid centers 
Vnikx) correspond to a nonzero Chern number. These 
Chern numbers can only be removed by breaking inver¬ 
sion symmetry of individual I'l/'nk)) similar to the TR 
symmetric case discussed in Sec. HB. Therefore, the ob¬ 
struction to smoothness persists in the inversion sym¬ 
metric gauge. Note that also in a Z 2 -even case inversion 
symmetry can protect additional topologies as has been 
pointed out in Ref. 36 and 46. 

Thus, the adiabatic connection of an inversion- 
symmetric topological insulator to a trivial one should 
be found by breaking the inversion symmetry along the 
path. 


IV. APPLICATION TO KANE-MELE MODEL 

We first illustrate our technique by applying it to the 
Kane-Mele (KM) model that describes a 2D quantum 
spin Hall (Z 2 -odd) insulator in some of its parameter 
space.*® A smooth gauge and the corresponding WFs 
were obtained previously for this model by other meth¬ 
ods,and it is instructive to validate our method ver¬ 
sus known results before going to more complicated cases. 
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FIG. 3. Phase diagram of the KM model for Aso = 0.6 and 
Ar = 0.5. The red arrows indicate the path used to construct 
an obstruction-free gapped path between the trivial and topo¬ 
logical phases. The symmetry breaking is indicated in the 
third, out of plane, axis. The BZ of a honeycomb lattice is 
shown in the inset. 

A. Kane-Mele model 

The KM is a tight-binding model on a honeycomb 
lattice with one spinor orbital per site. The primi¬ 
tive hexagonal lattice vectors are Oi = ax and 02 = 
^(x + ^y) with the atoms A and B located at at sites 
tyi = |(a.i + 0 , 2 ) and ts = |(£ii + 0 , 2 )- The KM Hamil¬ 
tonian is given by 

^ =Y1 + iAso 

T f Ar 'y ^ ^ T 

<ij> i 

where summation over the suppressed spin indices is as¬ 
sumed. The summation < ij > runs over all nearest 
neighbors, and the sum over <C ij runs over all second- 
nearest neighbors. Vij = (2/-\/3)(<ii x ^ 2 ) = ±1, with di 
and d 2 being the first-neighbor bond vectors encountered 
by an electron hopping from j to f, s is a spin-1/2 oper¬ 
ator. Aso and Ar are parameters defining the spin-orbit 
coupling and is a staggered onsite potential. In what 
follows we fix t = 1, Aso = 0.6 and Ar = 0.5. 

B. Constructing an obstruction-free path for the 
Kane-Mele model 

Following our method, first an adiabatic path connect¬ 
ing the normal insulator (NI) to the TI phase of the KM 
model needs to be found. Let us first tune Xi, while keep¬ 
ing the spin-orbit coupling parameters Aso and Ar fixed. 
This path corresponds to the vertical axis in Fig. 3. A 
gap closure occurs at the two Dirac points at K and K' in 
the BZ (illustrated in the inset of Fig. 3). We thus need 


to find a symmetry breaking held that would prevent the 
gap closure at both Dirac points simultaneously. 

We thus introduce a hopping anisotropy S to the 
nearest-neighbor hopping tij. For the reasons that are 
made clear below, it is chosen such that the hopping is 
t(l -b d) in the (1,1) direction and t in the other direc¬ 
tions, and it breaks the Ca-rotational symmetry of the 
KM model (see Appendix for a thorough discussion of 
the symmetries in the KM model). By tuning S and mov¬ 
ing along the A^ = 0 line in Fig. 3 the gap closure can be 
moved to the M point. This path is depicted by the red 
arrows in Fig. 3. 

The degeneracy at the M point is easily lifted by an ap¬ 
propriate weak TR-symmetry breaking field. We found 
the optimal field to be a staggered magnetic held in the 
(—1, l)-direction of the form®^ 

^{2) = + \^v + ( 14 ) 

where a and r are the Pauli matrices acting in the spin 
and sublattice subspaces correspondingly. Apart from 
the TR-symmetry, this held breaks the mirror and C 2 - 
rotational symmetries of the KM model (see Appendix ) 
and the resultant WFs coincide with the ones reported 
in Ref. 43. 

The hnal path, shown by the red line in Fig. 3, is thus 
made of two parts: the line PI P2, connecting two points 
with the same Z 2 invariant, and the line P2 P3, where the 
Z 2 invariant changes from even to odd, and along which 
the TR symmetry is broken. Note that it is also possible 
to start the path directly at P2, which corresponds to 
the NI phase. However, for reasons made clear below, we 
find it illustrative to add P1P2 to the discussion. 


C. Wannier functions of the KM model 

Now that the gapped path connecting points PI and 
P3 is found, the above outlined method is implemented 
and the WFs are calculated along this path. The position 
space density®^ of the WFs is shown in Fig. 4 correspond¬ 
ing to the points P1,P2 and P3 on the adiabatic path of 
Fig. 3. 

At PI (NI phase) both WFs are localized on the lower 
energy sites. Each of them is mapped onto the other by 
TR, forming a Kramers pair of WFs with opposite spins. 

The sum of the Wannier centers P = fi -b f 2 , be¬ 
ing proportional to the electronic polarization,^^ is illus¬ 
trated with a red dot in the left panel in Fig. 4. The 
Ca-symmetry of the KM model (see Appendix ) con¬ 
strains the possible values of this sum, to the values 
that are invariant under this symmetry modulo a lat¬ 
tice vector. Consistent with Ca are the following val¬ 
ues: (P,,Py) = (0,0),(l/3,l/3),(2/3,2/3). A change 
between distinct values cannot occur continuously, unless 
the Ca symmetry is broken. The value of P correspond¬ 
ing to the point PI is Px = 1/3 and Py = 1/3. This 












V. APPLICATION TO BbSea 



FIG. 4. Charge density distribution corresponding to the 
WFs obtained at points PI, P2 and P3 of the Fig. 3. The 
left and right panels show the two WFs of the KM model 
at each of the points. The red point marks the sum of the 
corresponding Wannier centers (modulo a lattice vector). 


is yet another example of additional topological protec¬ 
tion. If there is a symmetry that quantizes the values of 
electronic polarization^^, and if the values in the topolog¬ 
ically trivial and nontrivial phases are different, then the 
symmetry must be broken so that the polarization can be 
changed continuously between the two discrete values.®^ 

When going from PI to P2, d 0 and A,, 7 ^ 0 and both 
the C 3 and C 2 -symmetries are broken. This allows to 
continuously change the electronic polarization without 
closing the band gap. At P2 (NI phase) the two lattice 
sites have equal energies, since = 0. Accordingly, 
the two WFs are equally distributed between both sites, 
again forming a Kramers pair. Note that Px = Py = 0 
at this point. 

At P3 we are in the desired Z 2 -odd phase. The field 
of Eq. 14 is introduced to get from P2 to P3 along the 
gapped path. The two WFs are localized on different 
sites, having opposite spins, in accord with the earlier 
study of Ref. 43. Clearly, this configuration does not 
conserve the TR symmetry of the system. The sum of 
the Wannier centers remains at Px = Py = 0. 


The standard example®’®^ of a TI is Bi 2 Se 3 . It has a 
rhombohedral lattice structure with the space group 
The material is layered with hexagonal quintuple layers, 
consisting of three Se (two equivalent Sel and one Se2) 
and two equivalent Bi atoms (shown as A, B and C in 
Fig. 6 ) bounded together by van der Waals interaction. 
The structure is inversion-symmetric, with an inversion 
center at the central SeO atom marked with the cross in 
Fig. 6 . 

Without SOC the band structure of Bi 2 Se 3 is topo¬ 
logically trivial, turning SOC on with a parameter Aso 
drives it into the TI phase. At Aso = 0 Bi 2 Se 3 is a di¬ 
rect small band-gap semiconductor (NI phase), whereas 
at the full experimental SOC strength Aso = 1 it is a 
TI. At the topological phase transition there is an inter¬ 
mediate semi-metallic state with a gap closure. This gap 
closure can be avoided by applying a suitable symmetry 
breaking field. 


A. Constructing a model Hamiltonian for Bi2Se3 

A fully self-consistent density functional theory cal¬ 
culation was carried out for Bi 2 Se 3 without spin-orbit 
coupling (SOC) using the Vienna ab initio simula¬ 
tion package (VASP)®®’®® with the projector augmented- 
wave method,®^ using generalized gradient approxima¬ 
tion of the Perdew-Burke-Ernzerhof for the exchange- 
correlation potential.®® The WannierQO^® package was 
then used to first disentangle an isolated group of bands 
and then to project the band structure onto the atomic 
p-orbitals of all Bi and Se atoms. No further iterative 
minimization of the WF spread was done. 

The resultant band structures, obtained from the pro¬ 
jected Hamiltonian and the ab initio calculations are il¬ 
lustrated in the left panel of Fig. 5. SOC is added after¬ 
wards to this Hamiltonian by adding a local Hamiltonian 
jjSOC basis of the atomic p-orbitals.®®’^® The case 

of full SOC is illustrated in the right panel of Fig. 5. The 
correct topological phase of the projected Hamiltonian 
was confirmed using the method of Ref. 38 

In the following we work with the projected Hamil¬ 
tonian, which allows for the easy control of SOC and a 
symmetry breaking field in the total Hamiltonian 

Htot Xso + a FT®®. (15) 

The term iF®® is introduced below. Note that both, 
and Ff®® are local. The two parameters Aso and a 
are found to be sufficient to construct an adiabatic path 
between the NI and TI phase in Bi 2 Se 3 . 
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FIG. 5. Projected band structure (Se (Bi) character shown in red (blue)) and the ab initio band structures (black dotted lines) 
of Bi 2 Se 3 . Left panel: no spin-orbit coupling. Right panel: spin orbit coupling included. 



A(Se1) 

B(Bi1) 

C(SeO) 

A(Bi2) 

B(Se2) 


FIG. 6. The layered structure of Bi 2 Se 3 . The cross denotes 
the inversion center at the SeO atoms. The black arrows indi¬ 
cate the staggered magnetic field on the Bi sites corresponding 
to iLff- The red arrows indicate an additional field corre¬ 
sponding to H^fse- 



B. Constructing an obstruction-free path for 
Bi2Se3 

The SOC strength Ago is tuned from Ago = 0 (NI 
phase) to Ago = 1 (TI phase with full experimental SOC 
strength) passing through a topological phase transi¬ 
tion at Ago ~ 0.47. At this phase transition the gap 
closes and a 3D Dirac cone is formed at the T point. 
To obtain a gapped adiabatic path connecting NI and 
TI phases, a suitable symmetry breaking held that 
gaps the Dirac cone needs to be found by considering 
the symmetry breaking requirements. Due to the dou¬ 
ble protection described in Sec. Ill, the held must break 
both TR and inversion symmetries. This requirement is 
satished by a staggered magnetic held, analogous to that 
of Eq. (14) used in the case of the KM model. The op¬ 
timal direction of the staggered held was found to be in 
the plane of hexagonal layers (see Fig. 6). 

Examples of two such staggered helds are illustrated 
in Fig. 6. The hrst one, H^, shown with black arrows, 
acts on Bi sites only. While both TR and inversion sym¬ 
metries are broken by this held, the combination of the 
two symmetries survives. As a consequence, the band 
structure remains doubly degenerate upon the inclusion 
of this held. 


FIG. 7. The bulk energy gap as a function of Aso- The 
symmetry breaking is tuned as sin(Aso7r). The field I7|f 
preserves the double degeneracy of bands. The gap is larger 
for the iffile field, which lifts this degeneracy. 

The other symmetry breaking held, breaks this 

compound symmetry as well by applying, in addition to 
, a held on all Se sites, as marked by the red arrows 
in Fig. 6. Application of this held lifts the double degen¬ 
eracy of the bands. Both helds prevent the gap closure 
and allow for the adiabatic connection between the NI 
and the TI. 

The band gap as a function of Ago is shown in Fig. 7. 
The parameter a of the TR symmetry breaking held is 
tuned as a(Aso) = sin(Ago7r). We interpolate the path 
between Ago = 0 and Ago = 1 using 9 intermediate 
equidistant steps. In what follows the BZ of the Hamil¬ 
tonian 15 is discretized into a 16x16x16 k-mesh. 


C. Wannier functions of Bi2Se3 

We hrst need to hnd WFs for the initial step of the 
path, that is for Ago = 0. Suitable trial states \Ti) for the 
topologically trivial band structure can be guessed from 
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Aso 

= 0 

Aso 

= 1 



Q. 


n 

rrSB 

^Bi 

11.01 

0.46 

11.04 

0.70 

^SB 

-^BiSe 

11.01 

0.46 

11.04 

0.73 


TABLE I. The spread of the resultant WFs after minimiza¬ 
tion. 

the occupations of atoms shown with color in Fig. 5(a). 
The occupied subspace consists mainly of the Se p-states. 
Therefore, the Se p-orbitals (both up and down spin) are 
chosen to be the initial trial states. 

Throughout the path we monitor the minimum of 
detS'(k) and the Wannier spreads. Apart from the first 
projection, where we start with local projections, the de¬ 
terminant always stays reasonably close to 1. For this 
reason, the minimization of the Wannier spreads done 
at each step of the path, after projecting onto the WFs 
obtained at the previous step, results in only small im¬ 
provement in localization. 

The Wannier spreads in the NI and TI phase are given 
in Tab. I. As expected we find them to be larger in 
the TI than in the NI phase.Both symmetry breaking 
fields lead to nearly identical WFs, although, the result¬ 
ing spread is found to be slightly larger when the field 
-^BiSe break the symmetry. The smoothness of 

the resulting gauge is visible in Fig. 8 (see Fig. 1), where 
we show the flow of the hybrid Wannier centers in the 
kx = 0 plane. All centers are smooth and periodic in the 
BZ. 

Now that a smooth gauge for Bi 2 Se 3 is obtained, the 
corresponding Bloch states can be used as a set of k- 
dependent trial states Ti(k). These states are useful to 
quickly obtain a smooth gauge and ELWFs on a differ¬ 
ent A:-mesh without repeating the above procedure. In 
the spirit of Wannier interpolation,^^ the new Ti(k) are 
obtained by Fourier transforming the ELWFs on the new 
k-mesh. 

The obtained ELWFs also allowed us to find the local 
trial states that work for Bi 2 Se 3 , that result in well lo¬ 
calized WFs for different k- meshes. They are listed in 
Tab. II, where the site location (see Fig. 6), the orbital 
character and the spin direction are provided for each of 
the trial states. Apart from the two exclusions, the states 
come in time-reversal pairs listed in the same row of the 
table. The trial states have mostly Se-character, apart 
from the exclusions, which are listed in rows 4 and 7 of 
the table, and reflect the band inversion that occurs be¬ 
tween the Px states of Se and Bi in Bi 2 Se 3 at the F-point, 
which can be seen in Fig. 5(b). 

As discussed above, Bi 2 Se 3 is an example of a double 
symmetry protection of the band topology, where both 
inversion and TR are broken in the smooth gauge. This 
breaking of the two symmetries is clearly seen in the 
flow of hybrid Wannier centers shown in Fig. 8. Fig¬ 
ure 8(a) shows hybrid Wannier centers constructed from 
a symmetry-preserving gauge. The topological obstruc¬ 
tion is evident. In contrast. Fig. 8(b) illustrates the 


ISeO.Px.Ti) 

|SeO,Pa:, li) 

|SeO,Pa,t^) 

|Se0,p„, 

|SeO,p,,t,) 

|SeO,pj,l;^) 

-^(|Sel,pj,,tj,) + |Bil,pi,,ti,)) 

|Sel,p,„,l,„) 

|Sel,p„,ti:) 

|Sel,p„,li.) 

|Sel.p,,t,,) 

ISel.p^.la,) 

|Se2,pj,,tj,) 

;^(|Se2,pj,,li,) + |Bi2,Pi.,li:)) 

|Se2,p„,ti.) 

|Se2,P;,,li:) 

|Se2.p,,t,,) 

|Se2,p,.l,,) 


TABLE 11. The 18 trial states for the topological material 
Bi 2 Se 3 . Left and right column show potential time reversed 
partners. 




FIG. 8. Flow of the hybrid Wannier centers in the = 
0 plane in BLSea. Panel (a): the gauge respects TR and 
inversion symmetries. Panel (b): the smooth gauge used to 
obtain the ELWFs. 


case of the smooth gauge obtained above. The smooth 
gauge Wannier centers are neither inversion, nor TR- 
symmetric, but are smooth and periodic in the BZ, cor¬ 
responding to zero individual Chern numbers. 


D. Evaluation of the Chern-Simons 
magnetoelectric polarizability 6 ^cs 

We now proceed to calculating the geometrical con¬ 
tribution to the magnetoelectric effect. The numerical 
evaluation of the Chern-Simons magnetoelectric polariz¬ 
ability ^cs introduced in Sec. I is a tedious task.^° While 
a direct simulation of a material’s response to electro¬ 
magnetic fields can potentially be used to evaluate this 
term, such a calculation requires the use of large super¬ 
cells, which makes it inefficient and computationally ex¬ 
pensive. Several methods to compute Ocs from the bulk 
wave functions exist,but only one of them 
(Ref. 30) was applied to materials within an ab initio 
framework. Our calculation is based on the formalism 
developed in the latter work, but the actual implementa¬ 
tion of the formulas is also done directly in position space, 
which significantly improves convergence with respect to 
the k-mesh, compared to the k-space implementations 
used before. 

To calculate 0cs the Bloch wave functions obtained 
from the bulk calculation are used. The Berry connection 
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matrix for these states is defined as 


d 

~ ('^mk I ^ | '^nk) 


(16) 


which is computed using the smooth gauge, obtained 
above/^ Then 0cs is evaluated by integrating the Chern- 
Simons 3-form over the entire Brillouin zone^’^^ 


(17) 

where the summation over band and Cartesian indices 
is assumed. While this integral is gauge invariant mod¬ 
ulo 27r, a smooth gauge is assumed in the derivation of 
Eq. 17.^’^"^ Thus a smooth gauge is required for a mean¬ 
ingful evaluation of the integral. This is when the above 
technique for finding a smooth gauge becomes important. 

Equation (17) can be rewritten in terms of matrix el¬ 
ements of position operators evaluated with WFs in po¬ 
sition space. Using 


Scs = - 


47r 




'BZ 


(18) 

R 


the result reported previously^® can be obtained 

^ R 

- \ y^(0/|r^|Rm) (Rm|rj |Pn) (Pn|rfc[OQ ). (19) 

RP ^ 


This equation still assumes a smooth gauge, since when 
the WFs are not exponentially localized, surface terms 
should be included to account for the slow decay of the 
Wannier matrix elements.®® 

We now use the ELWFs obtained for Bi 2 Se 3 above, 
to compute the 0cs-term in this material. While both 
formulas Eq. (17) and Eq. (19) agree in the limit of an 
infinitely dense k-mesh, the use of the position space for¬ 
mula of Eq. (19) results in faster convergence, compared 
to the k-space formulation of Eq. (17).™ 

We used the trial states of Tab. 11 to get localized WFs 
for a variety of k-meshes up to 80x80x80. The scaling 
in the TI phase is shown in Fig. 9. Because of the slow 
convergence in k-mesh density, an extrapolation to the 
infinitely dense mesh is required. A linear extrapolation 
in (Afc)^ is done using the last data points, getting in all 
cases very close to the expected value of 0cs = This 
extrapolation choice is dictated by the error in numerical 
evaluation of the Aj matrices and the Wannier matrix 
elements (Om|rj|Rn), which is of order 

Exponential convergence can be achieved by di¬ 
rectly evaluating (Om|rj|Rn) in position space with the 
ELWFs.This is, however, computationally more ex¬ 
pensive and was done only for relatively coarse k-meshes 
up to 25x25x25. The green triangles in Fig. 9 illustrate 
the results obtained from this position space approach. 
Even with this coarse k-mesh a reliable extrapolation to 



FIG. 9. Convergence of 6cs in the TI for varying densities 
of k-meshes. Ak is the nearest-neighbor spacing on the grid. 
Convergence to 6cs = rr can be reached for mesh densities of 
order 80x80x80 ((Afc)^ = 0.0068nm“^). The black circles 
(red crosses) correspond to the results obtained by evaluating 
the Berry connection in k-space without (with) maximal lo¬ 
calization of the WFs resultant from the projections of Tab. II. 
The green triangles correspond to the results obtained by eval¬ 
uating {Om|rj|Rn) in position space. A linear extrapolation 
to the infinitely dense k-mesh was done using the last data 
points. 



FIG. 10. 6cs calculated along the adiabatic path connecting 
NI and TI phases. The full (empty) circles correspond to the 
symmetry breaking field H^f (RBiSe)-^° 


the expected value Ocs = can be done within a few 
percent accuracy. 

The 0cs"term was also computed along the adiabatic 
path connecting NI to TI for the two different symmetry 
breaking fields (Rff and RBiSe)- The results are shown 
in Fig. 10. The 0cs-term is found to be larger for the 
RgP field, which does not break the product symmetry 
of TR and inversion. 
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VI. CONCLUSION 


In this paper we established a general procedure for 
constructing ELWFs and smooth Bloch states to describe 
a group of bands with nontrivial topology. This is done 
by connecting the topologically nontrivial Hamiltonian to 
some trivial one by a gapped adiabatic path. Our tech¬ 
nique works for all symmetry-protected topologies, pro¬ 
vided that the net Chern number of the bands is zero. It 
was illustrated for Z 2 TIs in two dimensions using the ex¬ 
ample of the Kane-Mele tight-binding model, and the real 
three-dimensional topological insulator material Bi 2 Se 3 . 
We also introduced the concept of double symmetry pro¬ 
tection of nontrivial topology, which describes (ubiqui¬ 
tous) situations when there is more than one symmetry 
that presents the obstruction for choosing smooth Bloch 
states. We expect that this technique can also be general¬ 
ized to the case of Chern insulators to obtain numerically 
a smooth lattice representation for all the occupied Bloch 
states, apart from one that carries the net Chern number 
of the system, and hence cannot be made smooth. 

Finally, we described how the proposed scheme of con¬ 
structing a smooth gauge allows for the calculation of 
the Chern-Simons 0cS"term that captures the geometric 
contribution to the orbital magnetoelectric response of 
materials. A detailed discussion of the numerical imple¬ 
mentation of this calculation was provided, showing how 
to efficiently implement our technique to materials, where 
9cs is not quantized. This technique is especially useful 
in the presence of band topologies that do not result in 
a quantized value of the 0cs-term. More generally, the 
numerical construction of the smooth gauge for lattice 
models presented here can have broader applications in 
evaluation of various (band) geometric effects that do not 
have an immediate gauge-invariant formulation. 
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Appendix: Symmetries of the Kane-Mele model 


Here the symmetries of the Kane-Mele model discussed 
in the main text are derived explicitly. Using the basis 
set (|A t), |H t), |A 4 ,), |H D) the Hamiltonian Eq. (13) 


symmetry 

representation R 

k± — y ... 

AI 2 —^ ■ • ■ 

(t(“> 


fc2 

ki 

(t1“> 

i Tx 

-k2 

-ki 

C 2 

iOzTx 

-ki 

-k2 

C 3 

"^0 

k2 — ki 

-ki 

Ce 

(#(To-bi^jr, 

k2 

k2 — ki 


TABLE III. Point symmetries (relative to the origin) of the 
KM model for A,/ = <5 = 0. Both momentum transformations 
and matrix representations are given. 


can be written in k-space 

H{k) = t ao[Ta;eti{k) -|-Tj^et 2 (k)] 

+Aso cr^TzesoO<-) 

+AR{cra;[Ta;eRi(k) -|-TyeR2(k)] (A.l) 

-b aylTyensO^) + ra:eR4(k)]} 

— Ay (ToT^, 


with cr and r matrices acting in spin and lattice sub¬ 
spaces; (To being the unit matrix, and cr^^y^z {jx,y,z) - the 
Pauli matrices. 

The k-dependent coefficients in the Hamiltonian are 


en(k) = (l + 5)/i(k) + / 2 (k) + / 3 (k), 

et 2 (k) = (1 + d)gi(k) -b g2{k) + S'3(k), 
eso(k) = 2 [- sin(/ci) -b sin(fc 2 ) -b sin(fci 


em = -^5i(k) - ^g2(k) -b 53(k), 

eR2 = -^/i(k) - ^/2(k) + /3(k), 

eR3 = ^/i(k) - ^/2(k), 

eR4 = ^51 (k) - ^52 (k). 


^ 2 )], 


(A.2) 


where 


/i(k) = cos(^-b ^), 5i(k) = sin(^-b'if). 


/ 2 (k) = cos(- 2 ^ -b ^), g 2 {k) = sin(- 2 ^ -b ^), 
fsik) = cos{!f - 2 ^), 53 (k) = sin(^ - 2 ^), 

with ki and ^2 being the reduced coordinates in terms 
of the reciprocal lattice vectors bi = ~ 

62 = —^ky. 

We first consider the case of Ay = d = 0, corresponding 
to the point P3 in Fig. 3. In this case the 2D space group 
(wallpaper group) of the model is p&m in lUC notation 
(*632 in orbifold notation),^® and the corresponding 2D 
point group is Dg. The specific matrix representation of 
all the point group symmetries for this case is provided 
in Tab. HI. Each of the symmetries acts according to 

H{S{h,k 2 )) = RsH{h,k 2 )Rl 


(A.3) 
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where S is the symmetry operation and Rs is its matrix 
representation. 

The consistency of the matrix representation of the 
Tab. Ill can be checked by noting that the C 2 symmetry 
is the product of the two mirrors and and 

the Ca symmetry is the Cq rotation applied twice, as 
expected. Other reflection symmetries can be obtained 
by appropriate combinations of these operations. The 
relations expected for spin-1/2 systems, namely = 
Cl = Cg = Cg = — 1 are satisfied. Without the Rashba 
term, that is for Ar, = 0, there would be an additional 
inversion symmetry (TqTx- 

The symmetries corresponding to the different points 
in the path of Fig. 3 are summarized in Tab. IV. When 
both Aj, ^ 0 and S ^ 0, which is the case on the line 
connecting PI and P2 (P1P2) in Fig. 3), the sum of 
the Wannier centers P is allowed to change continuously 
along the lines Py = ^Px + , where m is an integer. 


The TR symmetry in the model is given by layT^lC^ 
with /C being complex conjugation, and its action in k- 
space given by ki —>■ —ki and ^2 —k 2 - Two TR- 

breaking fields, and P[^^ were considered in this 

work 

Hfl) = {-^CTx + 

SB / 3 1 . 

^(2) — + ^0-y + -tTzjTx- 

Of them, commutes with both mirror symmetries 

The field is applied along the line P 2 P 3 , 
where, according to Tab. IV, both these mirrors are pre¬ 
served. Therefore, the resulting WFs reflect this sym¬ 
metry and the minimization of their spread gets stuck in 
a local minimum. By adding the component, lead¬ 
ing to both mirrors are broken and we obtain the 

maximally localized WFs also found in Ref. 43. 
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